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ABS ALGORITHMS FOR LINEAR 
EQUATIONS AND ABSPACK 

E. SPEDICATO(*) - Z. XIA(**) - E. BODON(***) - A. DEL POPOLO( t ) 



Abstract - We present the main results obtained during a research on ABS methods 
in the framework of the project Analisi Numerica e Matematica Computazionale. 



1. The ABS algorithms 

ABS algorithms were introduced by Abaffy, Broyden and Spedicato (1984), 
to solve linear equations first in the form of the basic ABS class, later generalized 
as the scaled ABS class and applied to linear least squares, nonlinear equations 
and optimization problems, see e.g. the monographs by Abaffy and Spedicato 
(1989) and Zhang, Xia and Feng (1999), or the bibliography by Nicolai and 
Spedicato (1997) listing over 300 ABS papers. In this paper we consider some 
new results obtained in the framework of project Analisi Numerica e Matematica 
Computazionale, including the performance of several codes of ABSPACK, a 
FORTRAN package under development. 

For later reference, we recall the scaled ABS algorithm for solving the follow- 
ing determined or underdetermined linear system, where rank(A) is arbitrary 
and A T = (ai, . . . , a m ) 

Ax = b x e R n , b G R m , m<n (1) 

or 

ajx — bi = 0, i = 1, . . . , m (2) 
The scaled ABS algorithm 
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(A) Give x\ G R n arbitrary, Hi <E R n ' n nonsingular arbitrary, v\ G R m arbi- 
trary nonzero. Set i = 1. 

(B) Compute the residual = Axi — b. If = stop (xi solves the problem.) 
Otherwise compute = HiA T Vi. If Si ^ 0, then go to (C). If s, = 
and r = vjri = 0, then set x i+i = Xi, H i+1 = Hi (the z-th equation is 
redundant) and go to (F). Otherwise stop (the system has no solution). 

(C) Compute the search vector pi by 

Vl = HTzi (3) 
where 2, € R n is arbitrary save for the condition 

vfAHfzi ± (4) 

(D) Update the estimate of the solution by 

x i+ i = Xi - ctiPi, cti = vfn/vfApi. (5) 

(E) Update the matrix Hi by 

H i+1 =Hi- H l A T v l wfH l /wfH l A T v l (6) 
where Wi £ R n is arbitrary save for the condition 

wfHiA T Vi ji 0. (7) 

(F) If i = m, stop (x m+ i solves the system). Otherwise give Vi+i E R m arbi- 

trary linearly independent from v\ , . . . , Vi . Increment i by one and go to 
(B). 

Matrices Hi, which are generalizations of (oblique) projection matrices, have 
been named Abaffians at the First International Conference on ABS methods 
(Luoyang, China, 1991). There are alternative formulations of the scaled ABS 
algorithms, e.g. using vectors instead of the square matrix Hi, with possible 
advantages in storage and number of operations. In a next section we will show 
how they can be used to generate infinite iterative methods. 

The choice of the parameters Hi, , Zi, Wi determines particular methods. 
The basic ABS class is obtained by taking Vi = ej, as the i-th unit vector in R m . 

We recall some properties of the scaled ABS class, assuming that A has full 
rank. 

• Define Vi = (v\, . . . , Vi) and Wi — (w\, . . . , Wi). Then H i+1 A T Vi = and 
HT +1 Wi = 0. 
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• The vectors HiA T Vi, Hfwi are zero if and only if m, Wi are respectively 
linearly dependent from a\, . . . , a,_i, w\, . . . , Wi-\. 

• Define Pi = [p\, . . . ,pi) and A.- L — (a\, . . . , ai). Then the implicit fac- 
torization V^AfPi = Li holds, where Li is nonsingular lower triangular. 
Hence, if m = n, one obtains a semiexplicit factorization of the inverse, 
with P = P n , V = V n , L = L n 

A' 1 = PL- 1 V T . (8) 

For several choices of V the matrix L is diagonal, hence formula (8) gives a 
fully explicit factorization of the inverse as a byproduct of the ABS solution 
of a linear system. 

• The general solution of system (1) can be written as follows, with q G R n 
arbitrary 

x = x m+ i + H^ +1 q (10) 

• The Abaffian can be written in term of the parameter matrices as 

H i+1 = h 1 - n x A T v i {yvjExA T v i )- x wJn x . (11) 

Letting V — V m , W — W m , one can show that the parameter matrices Hi, V, W 
are admissible (i.e. condition (7) is satisfied) iff the matrix Q = V T AH^W is 
strongly nonsingular (i.e. it is LU factorizablc). This condition can be satisfied 
by exchanges of the columns of V or W . If Q is strongly nonsingular and we take, 
as is done in all algorithms so far considered, z; L — Wi, then condition (4) is also 
satisfied. Analysis of the conditions under which Q is not strongly nonsingular 
leads, when dealing with Krylov space methods in their ABS formulation, to 
a characterization of the topology of the starting points that can produce a 
breakdown (either a division by zero or a vanishing search direction) and to 
several ways of curing it, including those considered in the literature. 

Two subclasses of the scaled ABS class and particular algorithms are now 
recalled. 

(a) The conjugate direction subclass. This class is obtained by setting Vi = pi. 
It is well defined under the condition (sufficient but not necessary) that 
A is symmetric and positive definite. It contains the ABS versions of 
the Choleski, the Hcstcncs-Sticfcl and the Lanczos algorithms. This class 
generates all possible algorithms whose search directions are A-conjugate. 
If x\ = 0, the vector Xi + ± minimizes the energy (yl-weighted Euclidean) 
norm of the error over Span(pi, . . . ,pi) and the solution is approached 
monotonically from below in the energy norm. 
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The orthogonally scaled subclass. This class is obtained by setting V{ = Api. 
It is well defined if A has full column rank and remains well defined even if 
m is greater than n. It contains the ABS formulation of the QR algorithm 
(the so called implicit QR algorithm) , the GMRES and the conjugate resid- 
ual algorithms. The scaling vectors are orthogonal and the search vectors 
are j4 T A-conjugate. If x\ — 0, the vector Xi+\ minimizes the Euclidean 
norm of the residual over Span{p\, . . . ,Pi) and the solution is monotoni- 
cally approached from below in the residual norm. It can be shown that the 
methods in this class can be applied to overdetermined systems of m > n 
equations, where in n steps they obtain the solution in the least squares 
sense. 

The optimally stable subclass. This class is obtained by setting Vi — A~ T pi, 
the inverse disappearing in the actual recursions. The search vectors in 
this class are orthogonal. If x\ — 0, then the vector x i+i is the vector of 
least Euclidean norm over Span(pi, . . . ,pi) and the solution is approached 
monotonically from below in the Euclidean norm. The methods of Gram- 
Schmidt and of Craig belong to this subclass. The methods in this class 
have minimum error growth in the approximation to the solution according 
to a criterion by Broyden. 

The Huang algorithm is obtained by the choices Hi = I , Zi = to, = a^, Vi = 
d. A mathematically equivalent, but numerically more stable, formulation 
is the so called modified Huang algorithm (pi = Hi(Hiai) and = Hi — 
Pipf /pfpi). Huang algorithm generates search vectors that are orthogonal 
and identical with those obtained by the Gram-Schmidt procedure applied 
to the rows of A. If x\ = 0, then Xi+\ is the solution with least Euclidean 
norm of the first i equations. The solution x + with least Euclidean norm 
of the whole system is approached monotonically and from below by the 
sequence Xi. 

The implicit LU algorithm is given by the choices H\ = I, Zi = Wi = 
Vi = ej. It is well defined iff A is regular (i.e. all principal submatrices 
are nonsingular). Otherwise column pivoting has to be performed (or, if 
m = n, equation pivoting). The Abaffian has the following structure, with 

Ki e R n -^ 

"0 



implying that the matrix Pj is unit upper triangular, so that the implicit 
factorization A = LP^ 1 is of the LU type, with units on the diagonal. 
The algorithm requires for m — n, n 3 /3 multiplications plus lower order 
terms, the same cost of classical LU factorization or Gaussian elimination. 
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Storage requirement for Ki requires at most n 2 /4 positions, i.e. half the 
storage needed by Gaussian elimination and a fourth that needed by the 
LU factorization algorithm (assuming that A is not overwritten). Hence 
the implicit LU algorithm has same arithmetic cost but uses less memory 
than the most efficient classical methods. 

(f) The implicit LX algorithm, see Spedicato, Xia and Zhang (1997), is defined 
by the choices Hi = I, Vi = e^, Z{ — Wi — e^, where ki is an integer, 
1 < ki < n, such that e^i^a^ ^ 0. By a general property of the ABS 
class for A with full rank there is at least one index ki such that eJ.Hiai ^ 
0. For stability reasons we select ki such that | eJ.HiOi | is maximized. 
This algorithm has the same overhead and memory requirement as the 
implicit LU algorithm, but docs not require pivoting. Its computational 
performance is also superior and generally better than the performance 
of the classical LU factorization algorithm with row pivoting, as available 
for instance in LAPACK or MATLAB, sec Mirnia (1996). Therefore this 
algorithm can be considered as the most efficient general purpose linear 
solver not of the Strassen type. The implicit LX algorithm has also very 
important applications in a reformulation of the simplex method for the LP 
problem, sec Zhang, Xia and Feng (1999), where it leads to a reduction 
of storage up to a factor 8 and of multiplications up to one order for 
problems where there is a small number of degrees of freedom, with respect 
to implementations based upon the classical LU factorization. 

2. Solution of linear Diophantine equations 

One of our main results has been the derivation of ABS methods for lin- 
ear Diophantine equations. The ABS algorithm determines if the Diophantine 
system has an integer solution, computes a particular solution and provides a 
representation of all integer solutions. It is a generalization of a method proposed 
by Egervary (1955) for the particular case of a homogeneous system. 

Let Z be the set of all integers and consider the Diophantine linear system 
of equations 

Ax = b, xeZ n , Ae Z mxn , b e Z m , m<n. (13) 

It is intriguing that while thousands of papers have been written concerning 
nonlinear, usually polynomial, Diophantine equations in few variables, the gen- 
eral linear system has attracted much less attention. The single linear equation 
in n variables was first solved by Bertrand and Betti (1850). Egervary was prob- 
ably the first author dealing with a system (albeit only the homogeneous one). 
Several methods for the nonhomogeneous system have recently been proposed 
based mainly on reduction to canonical forms. 
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We recall some results from number theory. Let a and b be integers. If there 
is integer 7 so that b — 7a then we say that a divides b and write a\b, otherwise 
we write a J[ b. If a\, . . . , a n are integers, not all being zero, then the greatest 
common divisor (gcd) of these numbers is the greatest positive integer S which 
divides all ai, i = 1, . . . , n and we write S — gcd(ai, . . . , a n ). We note that 8 > 1 
and that 5 can be written as an integer linear combination of the ai, i.e. S = z T a 
for some z G R n . One can show that 6 is the least positive integer for which the 
equation a\X\ + . . . + a n x n = S has an integer solution. Now S plays a main role 
in the following 

Fundamental Theorem of the Linear Diophantine Equation 

Let ai,...,a n and b be integer numbers. Then the Diophantine linear equation 
a\X\ + . . . + a n x n = b has integer solutions if and only if gcd{ai 1 . . . , a n )\ b. In 
such a case if n > 1 then there is an infinite number of integer solutions. 

In order to hnd the general integer solution of the Diophantine equation a\X\ + 
. . . + a n x n = 6, the main step is to solve aiXi + . . . + a n x n = <5, where 5 = 
gcd(a\, . . . ,a n ), for a special integer solution. There exist several algorithms 
for this problem. The basic step is the computation of S and z, often done 
using the algorithm of Rosser (1941), which avoids a too rapid growth of the 
intermediate integers, and which terminates in polynomial time, as shown by 
Schrijvcr (1986). The scaled ABS algorithm can be applied to Diophantine 
equations via a special choice of its parameters, originating from the following 
considerations and Theorems. 

Suppose Xi is an integer vector. Since Xi + i = x; t — ctiPi, then Xi + i is integer 
if ai and pi are integers. If vj Api\(vfri), then ai is an integer. If Hi and are 
respectively an integer matrix and an integer vector, then pi = Hfzi is also an 
integer vector. Assume Hi is an integer matrix. From (6), if vj AHfwi divides 
all the components of HiA T Vi, then Hi + \ is an integer matrix. 

Conditions for the existence of an integer solution and determination of all 
integer solutions of the Diophantine system are given in the following theorems, 
generalizing the Fundamental Theorem, see Esmaeili, Mahdavi-Amiri and Spedi- 
cato (1999), or Fodor (1999) for a different proof under somewhat less general 
conditions. 

Theorem 1 Let A be full rank and suppose that the Diophantine system (13) is 
integerly solvable. Consider the Abaffians generated by the scaled ABS algorithm 
with the parameter choices: Hi is unimodular (i.e. both Hi and H^ 1 are integer 
matrices); for i = 1, . . . ,m, Wi is such thatwj HiA T Vi = 5i, Si = gcd(HiA T Vi). 
Then the following properties are true: 

(a) the Abaffians generated by the algorithm are well-defined and are integer 
matrices 
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(b) if y is a special integer solution of the first i equations, then any integer 
solution x of such equations can be written as x = y + Hf +1 q for some 
integer vector q. 



Theorem 2 Let A be full rank and consider the sequence of matrices Hi gen- 
erated by the scaled ABS algorithm with parameter choices as in Theorem 1. 
Let X\ in the scaled ABS algorithm be an arbitrary integer vector and let Z{ be 
such that zjHiA T Vi = gcd(HiA T Vi). Then system (13) has integer solutions iff 
gcd(HiA T Vi) divides vjri for i = 1, . . . , m. 

We can now state the scaled ABS algorithm for Diophantine equations. 

The ABS Algorithm for Diophantine Linear Equations 

(1) Choose x\ £ Z n , arbitrary, H\ £ z nxn , arbitrary unimodular. Let i = 1. 

(2) Compute Tj = vjri and Sj = HiA T Vi. 

(3) If (si — and n = 0) then let x i+1 = x i} H i+ i = H i} r i+1 = n and go to 
step (5) (the ith equation is redundant). If (si = and ^ 0) then Stop 
(the ith equation and hence the system is incompatible). 

(4) {si ^ 0} Compute Si = gcd(si) and p, = Hj ' Zi, where Zi £ Z n is an 
arbitrary integer vector satisfying zfsi — Si. If Si An then Stop (the 
system is integerly inconsistent), else Compute on = Ti/Si, let Xi+\ = 

Xi — OLiPi and update Hi by H i+ \ = Hi % w T H % A i v . ' where Wi £ R n is an 

arbitrary integer vector satisfying wfsi — Si. 

(5) If i = m then Stop (x m+ i is a solution) else let i = i + 1 and go to step 
(2). 



It follows from Theorem 1 that if there exists a solution for the system (13), 
then x — x m+ i + H^ +1 q, with arbitrary q £ Z n , provides all solutions of (13). 

Egervary's algorithm for homogeneous Diophantine systems corresponds to 
the choices Hi = I, x x = and Wi — z i7 for all i. Egervary claimed, without 
proof, that any set of n — m linearly independent rows of H m+ i form an integer 
basis for the general solution of the system. We have shown by a counterexample 
that Egervary's claim is not true in general; we have also provided an analysis 
of conditions under which m rows in H m+1 can be eliminated. 
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3. The generalized implicit LU subclass 

The generalized implicit LU (GILU) subclass is denned by taking v% = z% = 
Wi = ei and Hi arbitrary nonsingular. The GILU subclass is well defined iff the 
matrix AHj E m is strongly nonsingular, where E m — (ei, e m ). 

The well definiteness condition involves the matrix AHf = (H\a\, . . . , Hia m ) T . 
H? can be interpreted as a right scaling or right conditioning operator on A, 
acting in the same way on the different rows of A. If A is full rank but not 
regular, the well definiteness condition can be satisfied by simply taking Hi as a 
suitable permutation matrix. By this choice LU factorization with column piv- 
oting is imbedded in the GILU subclass. It can also be shown that all sequences 
Xi generated by the basic ABS class can be obtained a suitable choice of H\ in 
the GILU subclass. 

The given parameter choices imply the following structure for the Abaffian, 

(14) 

where S n ~t,n £ i£™ _ *>™ is full rank. The total number of multiplications is no 
more than n 3 + 0(n 2 ) for m = n, a substantial saving over algorithms in the 
basic ABS class which may require 3n 3 multiplications, after the parameters are 
given. 

The GILU subclass is related to a representation of the ABS class given 
in Abaffy and Spedicato (1989) in terms of n — i vectors as a generalization 
of a method proposed by Sloboda (1978). In such a representation one takes 
Wi = Zi and Hi arbitrary, assuming that the feasible parameters z^s are known 
initially. Then the search vector can be written in the form pi = u\ where 
Uj = Hi T Zj, j = 1, . . . , n and the vectors it*- are updated, for i = 1, . . . , m, by 

uf 1 =u\- (afui/afuM- 

The relation between the GILU subclass and the representation in terms of 
n — i vectors is given by the following 

Theorem 3 Define the matrix Ui = (ui, . . . , u n ) by Uk = for k < i, Uk — u\ 
for k > i. Then Ui = Hj ', u\ = pi where Hi, pi are respectively the i-th Abaffian 
and search vector of the GILU subclass with initial Abaffian H[ = Z T Hi, Z = 
(zi, . . .,z„). 

Further analysis shows that the matrix Hi needs not be given explicitly at the 
initial step of the algorithm. The i-th row of Hi can be given just at the i-th step. 
It can be considered as a vector parameter (right scaling parameter), arbitrary 
save that the matrix Af(H^)i must be nonsingular, where {Hj)i is the matrix 
comprising the first i columns of Hj. Therefore this formulation shows that all 
right prcconditioners can be imbedded in the ABS class, right preconditioning 
being just equivalent to a change in the initial matrix Hi. 
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4. ABS Methods for KT Equations. 



The KT (Kuhn- Tucker) equations are the following special linear system, 
related to the optimality conditions when minimizing a quadratic function with 
Hessian G G R n ' n subject to the linear equality constraint Cp = c, C G 
R m > n , p,ge R n , c,zeR m 



G C 1 
C 



(15) 



If G is nonsingular, the coefficient matrix is nonsingular iff CG~ X C T is non- 
singular. Usually G is nonsingular, symmetric and positive definite, but this 
assumption, required by several classical solvers, is not necessary for the ABS 
solvers. 

To derive ABS methods using the structure of system (15), observe that (15) 
is equivalent to the two subsystems 

Gp + C T z=g (16) 

Cp = c. (17) 
Consider the general solution of Cp = c in the ABS form, with q G R n arbitrary 

P = V m +i + Hj n+1 q (18) 

The parameters used to construct p m +i an d H m+ \ are arbitrary, hence (18) 
defines a class of algorithms. 

Since the KT equations have a unique solution, there is a q which makes 
p the unique n-dimensional subvector defined by the first n components of the 
solution of (15). By multiplying Gp + C T z = g on the left by H m+ i we obtain 
the equation 

H m+1 Gp = H m+1 g (19) 
which does not contain z. Now there are two possibilities for determining p: 

(Al) Consider the system formed by (19) and (17). Such a system is solvable 
but overdetermined. Since rank(i? m+ i) = n — m, m equations are recog- 
nized as dependent and are eliminated in step (B) of any ABS algorithm 
applied to this system, which then computes the unique solution. 

(A2) In equation (19) replace p by the general solution (18) to give 

H m +iGH^ +1 q = H rn+ ig - H m+ \Gp m+ i. (20) 

The above system can be solved by any ABS method for a particular solu- 
tion q, m equations being again removed at step (B) of the ABS algorithm 
as linearly dependent. 
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Once p is determined, one can determine z in two ways, namely: 

(Bl) Solve by any ABS method the overdetermined compatible system 

C T z = g-G P (21) 

by removing at step (B) of the ABS algorithm the n — m dependent equa- 
tions. 

(B2) Let P = (pi, . . . ,p m ) be the matrix whose columns are the search vectors 
generated on the system Cp = c. Now CP — i, with L nonsingular 
lower diagonal. Multiplying equation (21) on the left by P T we obtain a 
triangular system, defining z uniquely 

L T z = P T g - P T Gp. (22) 

Extensive numerical testing has evaluated the accuracy of the above ABS algo- 
rithms for KT equations for certain choices of the ABS parameters (correspond- 
ing to the implicit LU algorithm with row pivoting and the modified Huang 
algorithm). The methods have been tested against the method of Aasen and 
methods using the LU and the QR factorization. The experiments have shown 
that some ABS methods are the most accurate, in both residual and solution 
error; moreover some ABS algorithms are cheaper in storage and in overhead, 
up to one order, especially for the case when m is close to n. In particular two 
methods based upon the implicit LU algorithm not only have turned out to be 
more accurate, especially in residual error, than the method of Aasen and the 
method using QR factorization via Houscldcr matrices, but are also cheaper in 
number of operations (the method of Aasen has a lower storage for small m but 
a higher storage for large m). 

In many interior point methods the main computational cost is to compute 
the solution for a sequence of KT problems where only G, which is diagonal, 
changes. In such a case the ABS methods, which initially work on the matrix 
C, which is unchanged, have an advantage, particularly when m is large, where 
the dominant cubic term decreases with m and disappears for m = n, so that 
the overhead is dominated by second order terms. Again numerical experiments 
show that some ABS methods are more accurate than the classical ones. 



5. — A class of ABS methods for matrix equations 



It is common, in particular in optimization, to find systems of matrix equa- 
tions 

A i »X = b i , i = l,...,m (23) 
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where operation • is denned by A • B —tr{A T B). We can write systems (23), 
with obvious definition of o, as 

AoX = b (24) 
where A has the following form, with A 1 G R n ' n for i = 1, . . . , m 

" A 1 



A = 



(25) 



Problem (24) is a linear system in the space of matrices and the associated 
projection operators arc matrices whose elements are matrices in R n ' n . This 
observation led us to consider a linear space denoted by (ij"> n ) n > n and study its 
linear algebra. The isomorphisms between R n ' n and R n , between (ij n > n ) n > n and 
R n allow to establish ABS algorithms to solve (24), to generalize the Huang 
algorithm and implicit LU algorithm and to define other special algorithms. 
Quasi-Newton matrices satisfying linear relations (for example, symmetry and 
sparsity) can be described in the considered matrix form and can be solved by 
the proposed matrix ABS algorithm. 

The ABS method for finding a solution X G R n ' n , of system (24) is as follows, 
the symbol * indicating transposition in the matrix space. 

The matrix ABS algorithm 

Step 1 Give X 1 G A"-"™'", W 1 G (R) n ' n , set k = 1. 

Step 2 Compute r k = A k • X k - b k and S k =H k oA k . 

Step 3 If S k + go to Step 4; if S k = and r = set X k+1 = X k , H k+1 = U k 
and go to Step 7 if k < m; otherwise stop. If S k = and t^O stop, the 
system (24) is incompatible. 

Step 4 Compute P k E R n > n by 

pk = ( H ky Q Z k ^ 

where Z k G R n ' n is arbitrary save that A k • P k ^ 
Step 5 Update the approximation of the solution by 

X k+1 =X k -a k P k , a k =r k /A k .P k (27) 
If k = m stop; X m+1 solves the system. 
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Step 6 Update the matrix H k by 

n k+i = H k _ H k A k ^ W fcj. . ^]/v{/ fe . (ft fe o ,4 fc ) (28) 

where W k £ R n ' n is arbitrary save that W k • (H k o A k ) ^ 
Step 7 Increment the index k by one and goto Step 2. 

Properties of the matrix ABS method generalize properties of the scaled ABS 
class, albeit proofs are not always obvious, see Spedicato, Xia and Zhang (1999). 
We just recall that the general solution of (24) can be expressed as follows, with 
W e R n ' n arbitrary 

X = X m+1 + (H m+1 )* o W (29) 

There are natural generalizations of the Huang and the implicit LU algorithms 
in the ABS class. Huang matrix algorithm allows to construct solutions to the 
quasi-Newton equation 

B'5 = r (30) 

where 6 = x' — x and r = F(x') — F(x) when solving the nonlinear system of 
equations F(x) = 0, r = V f(x') — V/(x) when minimizing the unconstrained 
function f(x). Equation (30) can be solved also under the additional condition 
that some elements of the solution take prescribed values, by which way we can 
introduce the conditions of sparsity, symmetry and positive definiteness. 



6. — A class of ABS derived iterative methods 

When n is large, both storage and number of operations may be too large for an 
implementation of ABS methods using explicitly the AbafRans. One can develop 
methods with lower storage and operations by working with formulations of ABS 
methods that use vectors and then using restart or truncation. This approach 
leads to loss of termination and generates iterative methods. Here we describe 
one such class, which is derived from the following formulation of ABS methods 
in terms of k vectors at the k— th step, due originally to Bodon, see Abaffy and 
Spedicato (1989), and where we take Zk, Wk to be multiple of each other. 

The Bodon-ABS vector algorithm 

Let x x e R n be arbitrary. Let H x £ R n ' n and V = (v 1 , ...,«„)<= R n ' n 
be arbitrary 

nonsingular. 
For k = 1 to n 

Tk = v\Axk - v^b 
If k> 1 then 
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Pk = H T Z k 

For j = 1 to k — 1 

pl +1 = pi - ( v J A p{/ v J A Pj) Pj 

End 

Pk =Pk 

Else 

p k = H? z\ 
Endif 

a k = T k /(vlAp k ) 
Xk+i =x k ~ a k Pk 

End 

The above algorithm leads to iterative methods either via restart or via trunca- 
tion. Here we consider only truncation. If to is the number of available storage 
vectors, then we keep only information from the last to iterations, i.e. we replace 
all iterations from j = 1 to k — 1 by iterations from j = k — mtofc — 1. Strategies 
where the kept vectors are not necessarily the last m vectors may of course be 
considered. The matrix Hi should require low storage, hence we take Hi = I. 
Parameters Vk should be linearly independent and parameters Zk feasible (no 
division by zero). 

Algorithm ABS(m): the truncated Bodon-ABS vector algorithm 

Let xi e R n be arbitrary. 
Give the integer to, 1 < m < n. 
Do k = 1, 2, ■ • •, until convergence 
T fe = v^Axk - v^b 
If k > 1 then 

t = max(l, k — to) 

Pk = z k 

For j = t to k — 1 

Pl +1 = Pk - ( v j A Pk/ v j A Pj) Pi 
End 

Pk=p\ 
Else 

Pk = zi, (zi e R n ) 
Endif 

ak = Tk/iv^Apk) 
x k+ i =x k - a k pk 
End do 
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If m — 1 one can show that the truncated implicit LU, Huang and implicit 
QR algorithm generate respectively the Gauss-Seidel, the Kaczmarz and the De 
la Garza methods. For a comparison of storage and operations requirements 
with other well-know iterative methods see Spedicato and Li (1999). 

One can prove that algorithm ABS(to) is well-defined if, letting V k = (v t , ■ ■ ■ , v k ) 
and Z k = (z t , • • • , z k ), where t — max(l, k — to), then V k ' AZ k is strongly non- 
singular for all k. 

Without loss of generality, we can define the scaling parameters by 

v k = A- T Y Pk (31) 

where Y is a symmetric, positive definite matrix. The following choices for Y 
define, in the original scaled ABS class, the three subclasses considered in section 
1 and require a low storage for large sparse systems: Y = I (the optimally 
scaled subclass); Y = A T A (the orthogonally scaled subclass); Y = A T > 0, for 
A = A T > 0, (the conjugate direction subclass). 

Algorithm ABS(to) with the parameter choice (31) has variational properties, 
related to those of the original ABS class, namely that x k+ i minimizes the error 
F-weighted Euclidean norm over a linear variety spanned by the last to search 
vectors. 

One can also show that the truncated generalized conjugate direction al- 
gorithm of Dennis and Turner (1987) can be obtained by special choices of 
the parameters in Algorithm ABS(to). From this equivalence it follows that 
ORTHOMIN(m)), ORTHODIR(to)) and GMRES(m) are special cases of ABS(m). 

The convergence of ABS (to) at a linear rate can be proved by requiring that 
the angle between p k and the gradient r k of the Y— weighted error norm be 
uniformly bounded away from 90°. 



Theorem 4 . Suppose that there exists 7 > such that for the search vectors 
generated by Algorithm ABS(m) with choice (31) one has 

\Pkf k |> 7 II Pk h\\f k || 2 (32) 

for all k. Then the sequence x k converges to x* and satisfies 

2 

II x k+1 - x* \\ Y < (1 - C J J 1 / 2 || x k - X* \\ Y , (33) 

where Cond(Y) is defined by the ratio of the largest to the smallest eigenvalues 
ofY. 
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7. ABSPACK and its numerical performance 

The ABSPACK project aims at producing a mathematical package for solv- 
ing linear and nonlinear systems and optimization problems using the ABS al- 
gorithms. The project will take several years for completion, in view of the 
substantial work needed to test the alternative ways ABS methods can be im- 
plemented (via different linear algebra formulations of the process, different pos- 
sibilities of reprojections, different possible block formulations etc.) and of the 
necessity of comparing the produced software with the established packages in 
the market (e.g. MATLAB, LINPACK, LAPACK, UFO ...). It is expected that 
the software will be documented in a forthcoming monograph and willl be made 
available to general users. 

At the present state of the work FORTRAN 77 implementations have been 
made of several versions of the following ABS algorithms for solving linear sys- 
tems: 

1. The Huang and the modified Huang algorithms in two different linear 
algebra versions of the process, for solving determined, underdetermincd 
and overdetermined systems, for a solution of least Euclidean norm 

2. The implicit LU and implicit LX algorithm for determined, undcrdeter- 
mined and overdetermined linear systems, for a solution of basic type 

3. The implicit QR algorithm for determined, underdetermined and overde- 
termined linear systems, for a solution of basic type 

4. The above algorithms for some structured problems (including banded and 
block angular matrices) 

5. An ABS algorithm proposed by Adib and Mahdavi-Amiri (1999) equivalent 
to a block-two ABS method 

6. ABS versions of the GMRES method, some requiring less storage 

7. Several ABS methods for KT equations. 

For a full presentation of the above methods and their comparison with NAG, 
LAPACK, LINPACK and UFO codes see Bodon, Luksan and Spedicato (2000), 
Bodon and Spedicato (2000a, b). Some results are presented in the Appendix. 
There the columns refer respectively to: the problem, the dimension, the algo- 
rithm, the relative solution error (in Euclidean norm), the relative residual error 
in Euclidean norm (i.e. ratio of residual norm over norm of right hand side), the 
computed rank and the time in seconds. Computations have been performed in 
double precision on a Digital Alpha workstation with machine zero about 10~ 17 . 
All test problems have been generated with integer entries or powers of two 
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such that all entries are exactly represented in the machine and the right hand 
side can be computed exactly, so that the given solution is an exact solution 
of the problem as it is represented in the machine. Comparison is given with 
some LAPACK and LINPACK codes, including those based upon singular value 
decomposition (svd) and rank revealing QR factorization (gqr). 
Analysis of all obtained results indicates: 

1. Modified Huang is generally the most accurate ABS algorithm and com- 
pares in accuracy with the best LAPACK solvers based upon singular value 
factorization and rank revealing QR factorization; also the estimated ranks 
are usually the same. 

2. On problems whose numerical estimated rank is much less than the di- 
mension, one of the versions of modified Huang is much faster than the 
LAPACK codes using SVD or rank revealing QR factorization, even more 
than a factor 100. This is due to the fact that once an equation is recog- 
nized as dependent it docs not contribute to the general overhead in ABS 
algorithms. 

3. Modified Huang is generally faster and more accurate than other ABS 
methods and classical methods on KT equations. 

It should be noted that the performance of the considered ABS algorithms in 
term of times could be improved by developing block versions, as it is the case 
for the LAPACK codes, a work presently in progress. 
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8. Final remarks 

Additional work done in the framework of the project, not described here in 
detail, has involved the following topics. 

• Improvement of the performance of Newton method via a special trunca- 
tion approach, see Deng and Wang (1998). 

• Critical review of variable metric methods for unconstrained optimization 
with discussion of the ABS applications in this field, see Luksan and Spedi- 
cato (1998). 

• Analysis of the relations between the ABS methods and the classical method 
of averaging functional corrections, see Gredzhuk and Pctrina (1998). 

• Development of indefinitely preconditioned truncated Newton methods for 
large sparse equality constrained nonlinear programming problems, see 
Luksan and Vlcek (1998). 

• Computation and update of inertias of KKT matrices for use in quadratic 
programming, see Zhang (1999). 

• Further applications of the implicit LX algorithm to the simplex method 
for the LP problem, see Spedicato and Xia (1999). 

The field of ABS methods is now mature from a theoretical point of view, albeit 
there are exciting possibilities for applications to new fields, e.g. the eigenvalue 
problems. We expect that the completion of the project ABSPACK will provide 
a useful new instrument for users of mathematical software. 
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APPENDIX 

RESULTS ON DETERMINED LINEAR SYSTEMS 
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RESULTS ON OVERDETERMINED SYSTEMS 



Condition number: 0.16D+21 
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Condition number: 0.29D+18 
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RESULTS ON KT SYSTEMS 



Condition number: 0.26D+21 
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